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Abstract 

A  review  of  mathematical  models  of  lithium  and  nickel  battery  systems  developed  at  the  University  of  South  Carolina  is  presented.  Models 
of  Li/Li-ion  batteries  are  reviewed  that  simulated  the  behavior  of  single  electrode  particles,  single  electrodes,  full  cells  and  batteries  (sets  of 
full  cells)  under  a  variety  of  operating  conditions  (e.g.  constant  current  discharge,  pulse  discharge,  impedance  and  cyclic  voltammetry). 
Models  of  nickel  battery  systems  are  reviewed  that  simulate  the  performance  of  full  cells,  as  well  as  the  behavior  of  the  nickel  hydroxide  active 
material.  The  ability  of  these  models  to  predict  reality  is  demonstrated  by  frequent  comparisons  with  experimental  data. 

©  2002  Elsevier  Science  B.V.  All  rights  reserved. 

Keywords:  Simulations;  Impedance;  Porous  electrode  theory;  Continuum  models 


1.  Li-ion  battery  modeling 

1.1.  Background 

Among  the  rechargeable  battery  systems  available  at 
present,  the  Li/Li-ion  systems  provide  the  highest  energy 
densities  [1].  This  makes  the  Li/Li-ion  systems  the  most 
attractive  candidates  for  electric  vehicle  and  space  applica¬ 
tions.  However,  owing  to  various  practical  problems  asso¬ 
ciated  with  Li/Li-ion  cells,  their  commercial  use  is  restricted 
primarily  to  small  coin  cells  and  18650  cylindrical  cells. 
In  an  effort  to  make  the  large  scale  use  of  these  systems 
workable  and  safe,  mathematical  modeling  has  proved  to 
be  an  efficient,  cost  effective  tool.  The  modeling  efforts 
of  researchers  in  the  field  of  Li/Li-ion  battery  systems  are 
principally  based  on  the  isothermal  electrochemical  model 
developed  by  Doyle  et  al.  [1]  for  galvanostatic  discharge 
of  Li/Li-ion  cells.  In  this  work,  the  authors  have  developed 
a  very  general  one-dimensional  (ID)  model  applicable 
to  almost  any  of  the  existing  Li/Li-ion  systems.  Their  model 
uses  porous  electrode  theory,  a  macro-homogeneous  app¬ 
roach  developed  by  Newman  [2],  to  describe  the  potential 
variations  in  the  solid  and  solution  phases.  The  material 
balance  in  the  solution  phase  is  described  using  the  con¬ 
centrated  solution  theory  and  in  the  solid  phase  using  a 
Fickian  diffusion  equation  in  spherical  coordinates.  Confin¬ 
ing  only  to  isothermal  conditions,  the  authors  validate  their 
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model  by  demonstrating  good  agreement  with  the  experi¬ 
mental  data  [3].  Later,  Pals  and  Newman  [4],  Song  and 
Evans  [5]  and  Wang  et  al.  [6]  extended  this  model  to  include 
an  energy  balance  in  order  to  predict  cell  temperature.  Pals 
and  Newman  [4]  used  an  average  heat  generation  method 
while  Song  and  Evans  [5]  and  Wang  et  al.  [6]  used  local  heat 
generation  terms.  An  electrochemical-thermal  model  using 
local  heat  generation  terms  has  been  shown  to  be  more 
accurate  [7],  and  includes  all  the  features  of  the  isothermal 
model  developed  by  Doyle  et  al.  [1]  and  the  additional 
energy  balance  included  by  Song  and  Evans  [5]  and  Wang 
et  al.  [6].  Therefore,  it  is  deemed  important  at  this  point  to 
explain  the  development  of  such  a  model  for  a  typical  Li-ion 
system,  as  explained  in  the  work  of  Gomadam  et  al.  [7],  This 
model,  either  wholly  or  partly,  has  been  used  by  researchers 
for  modeling  single  electrode  particles,  single  electrodes, 
full  cells  and  batteries  (sets  of  full  cells)  to  simulate  various 
operating  conditions  like  cyclic  voltammetry,  constant  cur¬ 
rent  discharge,  pulse  discharge  and  impedance.  The  pro¬ 
blems  encountered  with  Li-ion  systems  like  poor  rate 
capability,  thermal  runaway,  occurrence  of  undesired  side 
reactions  and  capacity  fade  have  been  addressed  using  the 
model. 

1.2.  System  description 

A  schematic  representation  of  a  typical  Li-ion  cell  is 
shown  in  Fig.  1.  The  cell  sandwich  consists  of  five  regions: 
the  negative  electrode  current  collector  made  of  copper,  the 
porous  composite  negative  insertion  electrode,  the  porous 
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Nomenclature 

a  specific  area  of  an  electrode  (in  "  ) 

A  first  parameter  in  Eq.  (46) 

A,  coefficients  in  Eq.  (41) 

B  second  parameter  in  Eq.  (46) 

b  dimensionless  parameter  as  defined  in  Eq.  (40) 

c  concentration  in  a  phase  (mol/m3) 

Cp  volume  averaged  specific  heat  capacity 

(J/(kg  K)) 

D  diffusion  coefficient  in  a  phase  (m2/s) 

E0  formal  potential  corresponding  to  Eq.  (46)  (V ) 
Ea  Arrhenius  activation  energy  (J/mol) 

F  Faraday’s  constant  (96  487  C/mol) 

/  mean  molar  activity  coefficient  of  the  electro¬ 

lyte 

h  external  heat  transfer  coefficient  (W/(m2  K)) 

i  local  current  density  in  a  phase  (A/m2) 

io  exchange  current  density  at  an  interface 

(A/m2) 

k  rate  constant  (A/(m2  (mol/m3)3/2)) 

J  local  volumetric  transfer  current  density  due  to 

charge  transfer  (A/m3) 

L  length  of  the  cell  (m) 

m  outward  normal  to  a  boundary  (m) 

M  molecular  mass  (kg/mol) 

P  generic  property  defined  in  Eq.  (4) 

q  local  volumetric  heat  generation  (W/m3) 

/f  rate  of  reaction  (mol/(m3  s)) 

R  universal  gas  constant  (8.314  J/(mol  K)) 

r  radial  coordinate  (m) 

S  host  site 

t  time  (s) 

t+  transference  number  of  Li+  ions  in  solution 

v  specific  volume  for  packing  of  spheres 

Veq  equilibrium  cell  voltage  (V) 

x  coordinate  along  the  cell  length  (m) 

y  coordinate  along  the  cell  height  (m) 

z  coordinate  along  the  cell  width  (m) 

Greek  symbols 
a  transfer  coefficient 

d  radius  of  an  electrode  particle  (m) 

AH  enthalpy  of  reaction  (J/mol) 

£  volume  fraction  of  a  phase 

i]  surface  overpotential  (V) 

0  local  state  of  charge 

k  conductivity  of  the  electrolyte  (S/m) 

kd  diffusional  conductivity  of  the  electrolyte 

(A/m) 

1  thermal  conductivity  (W/(m  K)) 

p  density  (kg/m3) 

a  conductivity  of  the  matrix  phase  (S/m) 

(j)  local  potential  with  respect  to  Li/Li+  (V) 

X  volume  fraction 


Subscripts 

to  the  left  of  an  interface 
+  to  the  right  of  an  interface 

1  matrix  phase 

2  solution  phase  (except  in  Eqs.  (32)  and  (33) 
where  it  corresponds  to  the  lithium  deposition 
side  reaction) 

a  anodic 

amb  ambient 

app  applied 

c  cathodic 

Cu,n  copper  current  collector-negative  electrode 

interface 

Dx  corresponding  to  solid  phase  diffusivity 

Z)2  corresponding  to  solution  phase  diffusivity 

film  of  film 

H+  of  H+ 

i0  corresponding  to  equilibrium  exchange  current 

density 

j  defined  when  used 

L  when  larger  particles  are  controlling 

Li  of  lithium 

Li2C03  of  lithium  carbonate 

NiOOH  of  NiOOH 

Ni(OH)2  of  Ni(OH)2 

n  negative  electrode 

n,s  negative  electrode-separator  interface 

p  positive  electrode 

p,Al  positive  electrode-aluminum  current  collector 

interface 

ref  reference  state 

S  when  smaller  particles  are  controlling 

s  separator 

s,p  separator-positive  electrode  interface 

a  corresponding  to  solid  phase  conductivity 

k  corresponding  to  solution  phase  conductivity 

Superscripts 

0  initial 

eff  effective 

max  maximum 

s  surface 

brug  Bruggeman  exponent 


separator,  the  porous  composite  positive  insertion  electrode 
and  the  positive  electrode  current  collector  made  of  alumi¬ 
num.  The  composite  electrodes  are  made  of  their  active 
material  particles  (e.g.  MCMB  2528  in  the  case  of  negative 
and  LiCo02  in  the  case  of  positive),  held  together  with  a 
PVDF  binder  and  a  suitable  filler  material  such  as  carbon 
black.  From  its  fully  charged  state,  when  the  open  circuit 
potential  is  4.2  V,  the  discharge  behavior  of  the  cell  is 
predicted.  At  the  beginning  of  discharge,  the  negative 
electrode  is  fully  lithiated  while  the  positive  electrode  is 
ready  to  accept  lithium  ions.  During  discharge,  the  lithium 
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Composite  negative  Composite  positive 

electrode  (Anode)  Separator  electrode  (Cathode) 


Cu  current 
collector 


Active  Material 
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Al  current 
collector 


Electrolyte 


Active  material 
(LiCo02) 


Li-S  >Lf  +«-  +S 


Li'+e  +S  >  Li-S 
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&  Electrolyte 


x=0 
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Fig.  1.  A  schematic  of  a  typical  Li-ion  cell  as  well  as  a  porous  Ni-MH  or  Ni-Cd  cell.  In  the  case  of  nickel  system  replace  LiCoC>2  with  NiOOH,  LiPF6  with 
KOH,  and  Li^Ce  with  porous  MH  or  Cd.  O  =  (000). 


ions  deintercalate  from  the  negative  electrode  particles  and 
enter  the  solution  phase,  while  in  the  positive  electrode 
region  lithium  ions  in  the  solution  phase  intercalate  into  the 
LiCo02  particles.  This  results  in  a  concentration  gradient, 
which  drives  lithium  ions  from  the  negative  electrode  to  the 
positive  electrode.  The  cell  voltage  decreases  during  dis¬ 
charge,  as  the  equilibrium  potentials  of  the  two  electrodes 
are  strong  functions  of  the  concentrations  of  lithium  on  the 
surface  of  the  electrode  particles.  The  cell  is  considered  to 
have  reached  the  end  of  discharge  when  its  voltage  drops  to 


3.0  V.  The  kinetic,  ohmic  and  mass  transfer  resistances 
result  in  heat  generation  inside  the  cell.  This  heat  genera¬ 
tion  could  be  different  at  different  points  inside  the  cell. 
Thus,  there  develops  a  temperature  distribution  in  all 
directions  inside  the  cell,  in  addition  to  distributions  in 
potentials  and  concentrations.  A  complete  mathematical 
representation  of  such  a  system  should  include  equations 
that  describe:  (1)  mass  transport  of  lithium  in  the  solid 
phases,  (2)  mass  transport  of  lithium  ions  in  the  solution 
phase,  (3)  charge  transport  in  the  solid  phases,  (4)  charge 
transport  in  the  solution  phase  and  (5)  energy  transport  in 
the  whole  cell. 

1.3.  Generalized  model  development 
1.3.1.  Assumptions 

Listed  below  are  some  assumptions  made  to  develop  the 
model,  the  justification  and  usage  of  which  in  existing 
literature  are  cited  alongside. 


(a)  The  electrolyte  (e.g.  1  M  LiPF6  in  a  2:1  mixture  of  EC 
and  EMC)  is  assumed  to  be  binary  with  only  Li+  as  the 
electroactive  species  [1,2]. 

(b)  Side  reactions  are  neglected. 

(c)  Butler- Volmer  type  kinetic  expressions,  of  the  form 
described  below,  are  assumed  to  describe  the  charge 
transfer  processes  occurring  across  both  the  electrode/ 
electrolyte  interfaces.  Thus,  the  divergence  of  the  local 
solution  phase  current  density  at  any  point  is  given  by 
[2] 


otherwise 

In  addition,  the  variation  of  current  with  solid  and 
solution  phase  concentrations  is  accounted  for  by  using 
a  concentration-dependent  exchange  current  density 
[1]- 

ioj  =  kj(c™  -  cs/-(c°/-(c2)^,  j  n.p  (2) 

(d)  The  active  material  in  each  electrode  is  taken  to  be 
made  of  spherical  particles  of  a  single  size,  which 
enables  one  to  represent  mass  transport  inside  the 
particle  with  a  spherical  diffusion  equation.  Further,  the 
specific  areas  of  the  electrodes  can  be  calculated  from 
the  equation 

aj  =  3  ^ ,  j  =  n,  p  (3) 

°j 

(e)  The  reversible  part  of  the  heat  generated  due  to  charge 
transfer  reactions  at  the  interfaces  is  neglected.  This 
heat  is  proportional  to  the  quantities  d(f>jle{/dT,  which 


Vi  2  =  J  = 


exp 


JF 

RT 


exp 


JF 

RT 


1j 
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are  not  known  for  the  LirC6/LiPF6/LivCo02  system  and 
are  expected  to  be  small  [8].  Therefore,  the  variations 
of  (h;  ref  with  temperature  are  not  considered,  i.e. 

d<t>j,J/dT  =  0- 

(f)  Arrhenius  type  expressions  are  assumed  to  describe  the 
temperature  variations  of  kinetic  and  transport  para¬ 
meters  [6]. 


P  = 


Plef  exp 


£a,p 

(  1 

-Ml 

R 

\7ref 

t)\ 

r  k),n:  k),p:  0\n,  D i  p,  Z)2, 


(4) 


Other  parameters  like  f,  p,  Cp,  X  and  a  are  taken  to 
be  temperature  invariant. 

(g)  Volume  changes  arising  during  cell  operation  are 
neglected  resulting  in  constant  electrode  porosities  [1]. 

(h)  Constant  values  for  transference  number  and  solution 
phase  diffusivity  are  assumed  at  all  times  and  at  all 
points  in  the  cell. 


Expressions  obtained  from  fitting  the  experimental  data  are 
used  to  represent  the  variation  of  (f>jle f  with  9  [1].  The  last 
term  in  Eq.  (8)  is  included  in  order  to  correct  for  film 
resistances  at  the  electrode/electrolyte  interfaces. 

Applying  Ohm’s  law  for  charge  transport  in  the  solid  and 
solution  phases  after  modifying  it  to  account  for  concen¬ 
trated  solutions,  one  obtains  the  governing  equations  for 
potential  distributions  in  the  two  phases  [6]. 

V-(o-effV(/)1)-/  =  0  (11a) 


V  •  (fceffV(/>2)  +  V  •  (kdV lnc2)  +  J  =  0  (lib) 


where  the  effective  conductivities  are  given  by  Bruggeman’s 
correlation  [1] 


eff  brag,- 

K  =  1 , 


J  =  n,  s,p 


(12) 


and 


of 


=  tjjt: 


brug; 
IJ  > 


j  =  n,  P 


(13) 


1.3.2.  Governing  equations 

Application  of  energy  balance  [6,7]  around  an  arbitrary 
volume  element  inside  the  cell  yields 

dT 

pCp-=X7-(XX7T)  +  q  (5) 

where  the  local  source  term  is  given  by  [6]  (see  assumption 

(i)), 


q  =  trVc/q  •  S7(f)l  +  K'V (f>2  ■  V02  +  kd  •  Vlnc2  •  V02 

+  aJii(fj  +  /  =  n’p  (6) 

The  first  three  terms  arise  from  ohmic  heating  in  the  solid 
and  solution  phases.  The  last  term  is  the  heat  generated  due 
to  charge  transfer  at  the  electrode/electrolyte  interfaces.  This 
involves  a  reversible  part,  proportional  to  d(f>jlei/dT,  and  an 
irreversible  part,  proportional  to  rp.  Assumption  (e)  reduces 
the  heat  generation  to 

q  =  aX7(f)l  •  +  jcV<j)2  •  V02  +  KoVlnc2  •  V</)2 

+  ajij>1j,  ./  =  n,p  (7) 

where  the  surface  overpotential,  rjj,  is  defined  to  be 

>1j  =  01  -  02  -  ^'.ref  -  -%»  j  =  n>P  (8) 

Qj 

The  last  term  on  the  RHS  denotes  the  loss  due  to  a  resistive 
film  formed  over  the  electrode  particles.  The  value  of  the 
resistance  of  this  film  is  not  known  and,  therefore,  it  is  used 
as  an  adjustable  parameter.  The  equilibrium  potentials  (/>;,ref 
are  known  to  vary  strongly  with  state  of  charge  (SOC)  and 
are  expressed  as  functions  of  Qj  (  also  see  assumption  (f)) 
where 


9j  =  — — ,  /  =  n.  p 

J  ^max  ’  J  r 

ci  j 


(10) 


and  the  diffusional  conductivity  kd  is  given  by  [6] 


kd 


2  RTKeil(t+  -  1) 
F 


1  + 


din/' 

d  In  c2 


(14) 


From  assumption  (h),  the  differential  term  vanishes  from 
Eq.  ( 14).  Once  again,  the  experimental  data  has  been  fitted  to 
describe  the  variation  of  k  with  concentration  of  electrolyte. 

Species  conservation  on  any  electrode  particle  using 
assumption  (d)  results  in  the  Eq.  (1) 


dc  1  j  _  1  d  [  2  dcij 
dt  r2  dr  1,7  dr 


j  =  n,  p 


(15) 


and  on  the  solution  phase  results  in  the  equation  after 
assumption  (h) 

,2^=V.  (DfVc2)  +  ~t+)  J  (16) 

at  F 

The  effective  diffusivity  is  given  by  Bruggeman’s  correla¬ 
tion  [1]  as 

Df  =D22rj,  j  =  n,s,p  (17) 


1.3.3.  Initial,  boundary  and  interface  conditions 

Uniform  initial  conditions  were  used  for  T,  (p  and  c2 


at  t  = 

0, 

T  = 

rp 0 

for  all  x,y,z>  0 

(18) 

at  t  = 

0, 

CU~- 

for  all  x,y,z>  0 

(19) 

and 

at  t  = 

0, 

C2  = 

4 

for  all  x,  y,  z  >  0 

(20) 

At  all  boundaries,  except  the  current  collector/tab  inter¬ 
faces,  flux  boundary  conditions  are  applied  for  the  depen¬ 
dent  variables.  For  temperature,  the  flux  is  equated  to  the 
heat  lost  to  the  surroundings  according  to  Newton’s  law  of 
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cooling.  For  the  potentials  and  solution  phase  concentration, 
the  respective  fluxes  are  equated  to  zero.  Mathematically, 


dT 

-A—  =  h(T-T,mb) 
am 


(21) 


ftyi 

dm 


=  0 


(22) 


cells  [3].  Using  D2,  Dx,  Rf,  t+  and  brug  as  adjustable 
parameters,  they  could  obtain  good  agreement  for  low 
and  moderate  rates.  However,  in  a  later  communication, 
Arora  et  al.  [9]  showed  that  the  simulations  based  on  the 
same  model  do  not  agree  with  the  experimental  data  at 
higher  rates.  They  attributed  this  mismatch  to  one  or  more  of 
the  following  assumptions  involved  in  the  model: 


d<t>2  <91nc2 

K  ~r; - b  KD  — — —  0 

am  am 


(23) 


(24) 


where  m  denotes  the  outward  normal  to  the  boundary.  The 
boundary  condition  for  (f>i  changes  at  the  tab/current  col¬ 
lector  interfaces.  On  the  copper  current  collector/tab  inter¬ 
face,  <!>]  is  arbitrarily  set  to  0,  while  on  the  aluminum  current 
collector/tab  interface,  the  matrix  phase  current  density  is 
equated  to  the  applied  current  density,  i.e. 


d 

dm  ~  'app’ 


j  =  n,  p 


(25) 


For  diffusion  inside  the  electrode  particles  [1] 

A 

at  r  =  0,  =  0,  j  =  n,  p  (26) 


"  r  =  5”  =  J  =  n-P  <27) 

At  the  interfaces  between  the  different  regions  in  the  cell, 
the  following  conditions  are  applied  to  maintain  the  con¬ 
tinuity  of  fluxes.  At  all  interfaces,  all  fluxes  on  the  left  of  the 
interfaces  are  equated  to  those  on  the  right  [8].  The  excep¬ 
tions  are 


at  x  —  Lq L1  n, 

atx  =  Lns: 

at  x  =  Lvp . 


reff 

dc2 

dx 


d<t>  2 
dx 


.eff 


d  In  c2 


D  dx 


=  0 


=  0 


dh 

dx 

dh 

dx 


K 


at  x  =  L, 


=  0 

=  0 

b 

eff  d(f>2 

dx 


+< 


d  In  c2 


dx 


=  0 


-p,Al ; 


dc2 

dx 


=  0 


(28) 

(29) 

(30) 

(31) 


1.4.  Full  cell  isothermal  models 


Doyle  et  al.  [3]  validated  their  isothermal  model  by 
comparing  it  with  the  experimental  data  obtained  from 
Bellcore  (presently  Telcordia  Inc.)  plastic  Li-ion  (PLION) 


(i)  no  salt  precipitation; 

(ii)  no  solvent  transport  or  segregation  of  the  two  solvents; 
(iii)  constant  solid  and  solution  phase  transport  properties. 

They  investigate  the  reasons  one  after  another  in  view  of 
relaxing  them  to  obtain  better  fits  to  experimental  data. 

Salt  precipitation  may  occur  in  the  negative  electrode 
region  during  discharge  if  the  local  concentration  of  the 
electrolyte  reaches  values  above  its  solubility  limit.  For  the 
discharge  rates  considered  in  their  work,  they  showed  that 
the  electrolyte  never  reaches  saturation  and,  consequently, 
they  discarded  the  idea  of  including  salt  precipitation  effects 
in  their  model.  The  authors,  however,  deemed  solvent 
segregation  as  highly  possible  when  current  is  passed 
through  a  mixture  of  solvents.  Further,  the  polymer  electro¬ 
lyte  matrix  used  in  PLION  cells  was  treated  to  be  an  inert 
binder,  which  might,  in  reality,  interact  with  the  ionic 
species.  To  include  these  effects,  mass  transport  in  the 
solution  had  to  be  described  considering  the  electrolyte  as 
being  multi-component  in  nature  as  opposed  to  binary.  Such 
a  treatment  of  the  electrolyte  would  require  knowledge  of  10 
independent  transport  properties  for  the  five  species 
involved  (the  polymer,  the  two  solvents,  Li+  and  PF6“). 
Due  to  lack  of  available  data,  this  approach  was  also 
discarded. 

For  the  rates  (1.5-3C)  considered  in  their  work,  the 
authors  show  that  it  is  predominantly  the  solution  phase 
diffusion  that  limits  cell  performance.  Therefore,  it  was  clear 
that  the  value  of  D2  used  in  the  simulation  is  very  critical  in 
predicting  the  experimental  data.  Further,  the  fact  that  the 
local  solution  phase  concentration  varied  significantly  with 
position  inside  the  cell  at  high  rates  suggested  that  using  a 
concentration-dependent  D2  would  result  in  better  fits.  But 
since  no  experimental  data  was  available  relating  D2  with  c2, 
the  authors  tried  various  algebraic  expressions  of  the  form: 
(i)  D2  =  D2  +  ;«( 1  -  c2),  (ii)  Do  =  D2  exp (mc2),  (iii) 
Do  =  mf{K)  where  k  =f(c2),  etc.  They  used  each  of  the 
mentioned  expressions  in  the  simulations  but  found  that 
none  of  them  could  predict  the  experimental  data  any  better 
without  the  use  of  other  adjustable  parameters.  However,  by 
using  a  rate-dependent  D2  they  could  fit  all  discharge  rates 
well.  The  values  of  D2  used  for  three  different  cells  (called 
thin,  medium  and  thick  based  on  the  thickness  of  the 
electrodes  used)  are  shown  as  a  function  of  discharge  rate 
in  Fig.  2.  The  authors  acknowledge  that  this  makes  the 
model  empirical  and,  therefore,  usable  only  within  the  range 
of  discharge  rates  considered  in  their  work.  For  rates  as  high 
as  3C,  solid  phase  diffusion  limitations  were  shown  to 
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Fig.  2.  Solution  phase  diffusion  coefficient  as  a  function  of  discharge  rate  used  to  fit  experimental  data  for  three  different  cells.  1C  corresponds  to  1.156, 
1.937  and  2.691  A/m2  for  thin,  medium  and  thick  cells,  respectively. 

become  important  and  the  value  of  Dp  was  also  varied  with 
rate  to  obtain  good  fits  to  experimental  data. 

The  aforementioned  method  of  varying  D2  as  a  function  of 
rate  was  also  employed  by  Gomadam  et  al.  [10]  in  order  to 
obtain  good  fits  to  experimental  data  for  constant  current 
discharges  over  a  very  wide  range  of  rates  ((l/5)-10C). 

However,  unlike  Arora  et  al.  [9],  Gomadam  et  al.  [10]  found 
that  the  value  of  D2  had  to  be  increased  with  discharge  rate. 

Using  the  relation  between  D2  and  7ipp  from  constant  current 
discharges,  they  simulated  pulse  discharges  for  different 
amplitudes  (peak  currents),  duty  cycles  and  frequencies. 

The  simulated  cell  voltage  plotted  as  a  function  of  discharge 
capacity  in  Fig.  3,  compared  well  with  the  experimental 
data.  The  value  of  D2  during  the  intermittent  rest  periods  of  a 
pulse  discharge  was  used  as  an  adjustable  parameter  to 
obtain  good  fits  to  experimentally  obtained  data  for  voltage 
relaxation  of  the  cell  with  time.  It  should  be  noted  that  the 
authors  assumed  very  fast  transport  of  lithium  in  the  solid 
phase  at  all  rates. 

Guo  et  al.  [11]  developed  a  model  to  estimate  the  solid 
phase  diffusion  coefficient  from  experimental  impedance 
data.  Their  model  followed  the  model  of  Doyle  et  al.  [12]  for 
impedance  response  of  a  Li-ion  system.  Using  this  full  cell 
model  Guo  et  al.  [11]  simulated  the  impedance  response  of  a 
Li-ion  cell  as  shown  in  Fig.  4.  Also  shown  in  Fig.  4  are  the 
impedance  responses  of  the  individual  cell  components. 

Treating  these  as  experimental  data,  the  authors  described 
a  method,  called  the  modified  EIS  method,  to  extract  the 
solid  phase  diffusion  coefficient  by  the  technique  of  non¬ 


linear  parameter  estimation.  They  concluded  that  the  valid¬ 
ity  of  estimating  the  diffusion  coefficient  from  the  impe¬ 
dance  response  of  a  porous  intercalation  electrode  by  the 
modified  EIS  method  is  not  assured  if  the  transport  of 
lithium  ions  in  the  solution  phase  is  limiting  cell  perfor¬ 
mance.  They  identified  the  ratio  of  the  time  constants  of 
solid  and  solution  phase  diffusions  as  a  dimensionless 
parameter  governing  the  reliability  of  the  estimation  of 
the  solid  phase  diffusion  coefficient.  Their  results  showed 
that  bigger  particle  sizes,  smaller  electrode  thicknesses, 
larger  differences  in  the  values  of  the  solid  and  solution 
phase  diffusion  time  constants  and  small  active  material 
loading  are  more  conducive  to  a  valid  estimation  of  the  solid 
phase  diffusion  coefficient. 

The  model  developed  by  Doyle  et  al.  [1]  for  isothermal 
constant  current  discharge  of  Li-ion  cells  was  extended  by 
Arora  et  al.  [13]  to  include  the  lithium  deposition  side  reaction 
occurring  during  overcharge  conditions.  They  proposed  a 
simple  mechanism  of  Li+  +  e  — >  Li  for  the  reduction  of 
lithium  ions  on  parts  of  the  carbon  electrode  where  the 
potential  reaches  0  V  with  respect  to  a  Li/Li+  reference 
electrode.  The  total  local  reaction  rate  (7)  in  the  negative 
electrode  was  calculated  as  the  sum  of  two  individual  reaction 
rates:  (1)  the  rate  of  the  desired  lithium  intercalation  reaction 
(Jn  i)  and  (2)  the  rate  of  deposition  of  lithium  metal  (Jn  2). 

2 

jn  =  [exp(aa,k/>L)  -  exp(-ac,jl/t/A.)]  (32) 
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Fig.  3.  Comparison  of  simulated  voltage  profiles  with  experimental  data  for  pulse  discharge.  The  solid  lines  denote  model  predictions  while  the  circles 
denote  the  experimental  data. 
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where  «0,i  is  defined  by  Eq.  (2)  and 


io,2  =  k2c 2a'2  (33) 

The  authors  study  the  effect  of  lithium  deposition  on  charge 
efficiency  and  capacity  fade.  Further,  they  use  the  model  to 
optimize  certain  design  parameters  (like  negative  electrode 
thickness,  particle  size  and  excess  capacity)  and  operating 
conditions  (like  charging  rate  and  charge  cut-off  voltage). 
They  found  that  a  cell  with  5%  excess  negative  electrode 
when  charged  to  4.25  V  at  1C  rate  will  lead  to  lithium 
deposition,  which  occurs  after  the  cell  crosses  4.2  V.  Further, 
if  the  same  cell  had  19%  excess  capacity  in  the  negative 
electrode,  no  lithium  deposition  would  be  observed  even  if 
the  cell  was  overcharged  to  4.45  V  at  1C  rate.  However,  at 
higher  rates,  a  small  amount  of  lithium  deposition  will  be 
observed  when  the  cell  is  overcharged  to  4.45  V.  Thus,  the 
authors  conclude  that  the  lithium  deposition  side  reaction 
can  be  reduced  significantly  by  using  excess  negative  elec¬ 
trode  capacity,  although  the  downside  being  the  resulting 
enhanced  first  cycle  capacity  loss. 

The  authors  also  found  that  active  material  particle  size 
and  electrode  thickness  also  affect  lithium  deposition.  Thin¬ 
ner  electrodes  with  smaller  particle  sizes  were  found  to  be 
less  prone  to  lithium  deposition  than  thicker  electrodes  with 
larger  particles.  The  lithium  metal  deposited  on  the  negative 
electrode  was  also  considered  to  react  irreversibly  with  the 
solvent  forming  a  film  of  insoluble  products  (e.g.  LHCOs) 
around  the  active  material  particles.  This  resulted  in  an 
increase  in  the  cell  resistance,  which  was  shown  to  affect 
the  cell  capacity  loss  depending  on  its  composition.  The 
resistance  offered  by  this  film  to  current  flow  was  calculated 
as 

+  Xu  (— )  +  Zli.co,  (—)  (34) 

V  Ku  J  V^LbC 'Oj 

where  the  hint  thickness  (5film  is  given  by 

Ocifilm  ■  Al  2  H 

dt  anpF 

The  resistance  offered  by  the  him  was  found  to  increase  with 
an  increase  in  Li2C03  composition  resulting  in  earlier  end  of 
discharge  and,  consequently,  loss  of  discharge  capacity. 

The  models  discussed  thus  far  and  those  discussed  in  the 
subsequent  sections,  assume  that  the  electrodes  consist  of 
active  material  particles  of  one  size.  However,  in  reality, 
electrodes  contain  particles  of  a  range  of  sizes.  In  order  to 
understand  the  effect  of  multiple-sized  particles  on  the 
performance  of  the  cell,  Nagarajan  and  Van  Zee  [14]  devel¬ 
oped  a  model  assuming  that  the  cathode  consisted  of  par¬ 
ticles  of  two  different  sizes.  They  used  the  dilute  solution 
theory  to  describe  species  transport  in  the  solution  phase. 
Based  on  the  linear  packaging  model  developed  by  Yu  et  al. 
[15],  they  accounted  for  the  porosity  of  the  electrode  as  a 
function  of  the  mixture  of  particles  of  different  sizes.  The 
authors  considered  a  binary  mixture  of  small  and  large 


particles,  where  the  volume  fractions  of  the  small  (ys) 
and  large  (xl)  particles  are  known.  The  specific  volume 
of  the  mixture  is  calculated  based  on  the  volume  fractions  of 
the  particles  and  the  specific  volume  for  a  random  packing  of 
spheres  (v).  The  equation  for  the  specific  volume  will  depend 
on  the  array  matrix 

vl  =  vxL  +  vsXs  (36) 

i;s  =  vXs  +  vlXl  (37) 

Eq.  (36)  applies  for  a  matrix  of  large  particles  where  the 
voids  are  filled  with  small  particles,  while  Eq.  (37)  applies 
for  a  matrix  of  small  particles  where  the  voids  are  filled  with 
large  particles.  The  specific  volume  of  the  mixture  of 
particles  is  given  by  the  maximum  between  vs  and  vL. 
The  maximum  was  chosen  to  avoid  overlapping  of  the 
particles.  The  other  equations  used  to  complete  their  model 
are  given  as  Eqs.  (23),  (24),  (26)  and  (27)  in  their  paper  [14]. 
They  concluded  that  a  binary  mixture  can  provide  a  sig¬ 
nificantly  higher  density  of  active  material  relative  to  an 
electrode  comprised  of  single  size  particles.  However, 
increasing  the  packing  density  increases  the  liquid  phase 
diffusion  resistance,  and  therefore,  there  must  be  an  opti¬ 
mum  size  ratio  for  these  variables  that  maximizes  the 
capacity  of  the  cell. 

7.5.  Full  cell  non-isothennci l  models 


Pals  and  Newman  [4]  included  an  energy  balance  in  the 
isothermal  model  developed  by  Doyle  et  al.  [1]  in  order  to 
describe  non-isothermal  conditions  in  single  Li-ion  cells  and 
batteries  (sets  of  single  cells).  However,  they  used  a  lumped 
heat  generation  term,  based  on  the  difference  between  the 
equilibrium  cell  voltage  and  the  actual  cell  voltage,  to 
calculate  the  temperature. 


Veq  -  V  -  T(dVeq/dT)  . 

cl  —  .  'app 

^cell 


(38) 


Elsewhere,  Rao  and  Newman  [16]  showed  that  this 
method  of  calculating  the  heat  generation  term  is  accurate 
only  under  conditions  of  uniform  reactions  distribution  (e.g. 
low  rate  discharges).  Extending  the  single  cell  thermal 
model  developed  by  Pals  and  Newman  [4],  Botte  et  al. 
[8]  included  a  chemical  side  reaction  (the  negative  electrode 
decomposition),  which  being  an  additional  source  of  heat 
resulted  in  enhanced  temperature  rise  in  the  cells. 


V  -  T(dVeq/dr) 
7-cel! 


app 


-  AH[xnR 


(39) 


In  the  absence  of  such  side  reactions,  the  cell  temperature 
would  increase  only  due  to  the  irreversible  heat  generated  by 
electrochemical  charge  and  mass  transport  and  once  the  cell 
temperature  reaches  the  melting  point  of  the  separator  the 
cell  would  shut  down.  No  further  heat  will  be  generated  and, 
therefore,  the  cell  temperature  would  not  rise  above  the 
melting  point  of  the  separator.  On  the  other  hand,  if  there 
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were  chemical  side  reactions  occurring  inside  the  cell,  the 
cell  temperature  would  continue  to  rise  beyond  the  melting 
point  of  the  separator.  The  final  temperature  attained  by  the 
cell  has  been  shown  by  Botte  et  al.  [8]  to  be  sensitive  to  the 
reaction  parameters  namely  Ea  and  A//,xn.  They  also  studied 
the  effect  of  some  design  and  operating  parameters  on  the 
temperature  attained  by  the  cell  in  the  absence  of  any  side 
reactions.  The  heat  transfer  coefficient  has  been  pointed  out 
to  be  a  very  important  parameter  that  determines  the  heat 
generated  in  the  cell.  Starting  from  the  ambient  temperature 
of  25  °C,  the  authors  show  that,  for  a  3C  discharge,  the  cell 
temperature  reaches  values  as  high  as  90  °C  when  li  is  0  W/ 
(m2  K)  as  opposed  to  45  °C  when  h  is  20  W/(m2  K).  They 
also  show  that  the  final  cell  temperature  increases  from  30 
to  80  °C  when  the  rate  of  discharge  is  increased  from  1/2  to 
3 C  with  a  heat  transfer  coefficient  of  5  W/(m2  K).  Further, 
when  the  capacity  of  the  cell  is  increased  by  increasing 
the  electrode  thicknesses,  the  final  temperature  attained  by 
the  cell  increases  for  a  given  C-rate  discharge.  They  also 
studied  the  effect  of  electrode  particle  size  on  the  tempera¬ 
ture  attained  by  the  cell.  They  found  that  a  cell  with  30  pm 
particles  reached  60  °C,  while  that  with  10  pm  particles 
reached  only  40  °C  during  the  1C  discharge. 

Following  Wang  et  al.  [6],  Gomadam  et  al.  [7]  modified 
the  electrochemical-thermal  model  developed  by  Pals  and 
Newman  [4]  by  calculating  heat  generation  locally  instead 
of  using  a  lumped  method  (see  Eq.  (6)).  They  established  the 
validity  of  this  method  by  comparing  (see  Fig.  5)  the  single 
cell  model  predictions  with  experimental  data  obtained  from 
cylindrical  18650  cells.  The  voltage-time  predictions  of  the 
model  agreed  well  with  experimental  data  both  quantita¬ 
tively  and  qualitatively.  The  temperature-time  predictions, 
however,  agreed  only  qualitatively  well  and  showed  some 
quantitative  mismatches.  The  authors  attributed  this  to  errors 
in  estimation  of  unknown  parameter  values  and/or  violation 
of  one  or  more  of  their  assumptions  (e.g.  entropies  are 
negligible).  The  authors  also  extend  the  model  to  two 
dimensions  in  order  to  capture  non-uniformities  along  the 
height  of  the  cell.  They  show  that  these  effects  are  important 
particularly  when  the  height  of  the  cell  is  much  larger  than 
its  thickness.  The  final  temperatures  attained  by  the  top  and 
bottom  of  the  cell  was  found  to  first  increase  with  discharge 
rate  and  then  decrease  with  further  increase  in  the  discharge 
rate.  This  is  explained  based  on  the  fact  that  the  temperature 
attained  by  the  cell  is  proportional  to  the  product  of  the  rate 
of  heat  generation  (q)  and  the  discharge  time  (tD).  Since  q 
increases  and  rD  decreases  with  increase  in  discharge  rate, 
their  product  goes  through  a  maximum  and  the  cell  tem¬ 
perature  follows  suit.  The  authors  also  conducted  some  case 
studies  that  produced  counter-intuitive  results,  which  could 
be  predicted  only  by  calculating  the  heat  generation  term 
locally.  A  cell  in  which  current  enters  and  leaves  only 
through  the  tabs  connected  to  the  top  of  the  current  collector 
was  considered  to  lose  heat  to  the  ambiance  only  through  its 
bottom  face,  characterized  by  a  heat  transfer  coefficient  (h). 
In  such  a  case,  one  would  expect  the  cell  temperature  to 


decrease  at  all  points  when  the  value  of  heat  transfer 
coefficient  is  increased.  The  model,  however,  predicts  that 
the  temperature  at  the  top  portion  of  the  cell  would  increase 
while  that  at  the  bottom  part  would  decrease  as  expected. 
This  effect  has  been  explained  based  on  a  non-uniform 
distribution  of  heat  generation  (q)  inside  the  cell.  For 
non-zero  values  of  h,  as  discharge  proceeds,  the  bottom 
part  of  the  cell  loses  heat  and  hence  becomes  more  resistive 
while  the  top  part  becomes  hotter  and  more  conductive.  This 
makes  the  current  distribution  and,  consequently,  the  heat 
generation  even  more  non-uniform.  More  current  flows 
through  the  top  and  thus  the  top  part  gets  progressively 
hotter.  The  increase  in  the  local  solution  phase  current  (z2) 
near  the  top  part  of  the  separator  with  time  and  the  corre¬ 
sponding  decrease  in  i2  near  the  bottom  part  are  shown  in 
Fig.  6  for  a  cell  419  pm  thick  and  1  m  tall.  The  authors 
emphasize  that  the  lumped  heat  generation  method  used  by 
Pals  and  Newman  [4]  or  the  method  used  by  Song  and  Evans 
[5]  cannot  predict  such  results.  Further,  they  extend  the  two- 
dimensional  (2D)  single  cell  model  to  a  battery  (a  set  of 
single  cells)  using  the  modified  heat  generation  method 
explained  by  Pals  and  Newman  [4].  The  advantage  of  this 
method  is  that  the  electrochemical-thermal  model  has  to  be 
solved  only  sequentially  as  opposed  to  simultaneously  in  the 
rigorous  way. 

In  the  models  discussed  earlier,  the  intercalation  of 
lithium  into  the  solid  electrode  particles  was  described  by 
a  spherical  diffusion  equation  characterized  by  a  constant 
diffusion  coefficient.  However,  it  has  been  shown  that  the 
diffusion  of  lithium  inside  the  electrode  particles  is  non¬ 
ideal  and  does  not  obey  Fick’s  law  of  diffusion.  In  order  to 
accurately  describe  the  transport  of  lithium  through  the 
solid,  one  has  to  correct  the  diffusion  coefficient  by  an 
activity  term.  Verbrugge  and  Koch  [17]  developed  a  theo¬ 
retical  approach  to  calculate  this  activity  term  for  carbon 
fiber  electrodes  from  experimentally  measured  equilibrium 
potential  data.  Following  them,  Botte  and  White  [18] 
showed  the  importance  of  using  such  an  activity  corrected 
diffusion  coefficient  over  the  constant  diffusion  coefficient 
when  simulating  constant  current  discharges.  For  this  study, 
they  considered  a  cell  made  of  a  lithium  metal  cathode  and  a 
porous  carbon  fiber  anode,  which  are  the  same  as  that 
considered  by  Verbrugge  and  Koch  [17].  Following  the 
model  developed  by  Doyle  et  al.  [1]  for  constant  current 
discharges,  Botte  and  White  [18]  showed  that  for  low 
discharge  rates,  when  solid  phase  diffusion,  as  well  as  other 
transport  processes  do  not  limit  the  cell  performance,  the 
two  methods  predicted  identical  results.  However,  for  high 
rates,  the  constant  diffusion  coefficient  method  could  predict 
voltages  up  to  0.25  V  lower  than  those  predicted  using 
the  corrected  diffusion  coefficient  method.  Further,  the 
constant  diffusion  method  predicted  an  earlier  end  of  dis¬ 
charge.  The  authors  also  compared  the  cell  temperatures  as  a 
function  of  time  and  discharge  rate.  At  low  rates,  the  two 
methods  predicted  exactly  the  same  temperature  profiles 
while  at  higher  rates,  the  cell  temperature  as  predicted  by  the 
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Fig.  5.  Comparison  between  (a)  the  V-t  and  (b)  the  T—t  curves  obtained  for  different  discharge  rates  from  experiment  and  those  predicted  by  the  model.  The 
solid  lines  denote  model  predictions  while  the  circles  denote  the  experimental  data  (reproduced  from  [7]  by  permission  of  The  Electrochemical  Society). 


constant  diffusion  method  increased  more  steeply  with  time. 
However,  the  final  temperature  attained  by  the  cell  was 
higher  in  the  case  of  the  corrected  diffusion  method  since 
discharge  lasted  longer.  They  also  report  significant  mis¬ 
matches  between  the  two  methods  with  regard  to  the  con¬ 
centration  profiles  inside  the  solid  particles.  As  a  result  of  the 
aforementioned  observations,  the  authors  conclude  that  the 
corrected  diffusion  coefficient  method  must  be  preferred 


over  the  constant  diffusion  coefficient  method  especially  at 
high  discharge  rates. 

1.6.  Single  particle  models 

Zhang  et  al.  [19]  developed  a  spherical  diffusion  model  to 
predict  the  potentiodynamic  response  of  a  single  spinel  par¬ 
ticle.  Assuming  a  constant  solid  phase  diffusion  coefficient 
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Fig.  6.  Time  evolution  of  the  local  solution  phase  current  density  in  the  separator  at  the  top  and  bottom  of  the  cell,  h  —  104  W/(m2  K),  L  —  419  pm, 
/app  =  125  A/m2,  H  =  1  m,  T°  =  298  K,  ramb  =  0  K. 


(Dp)  and  using  that  as  an  adjustable  parameter  along  with  the 
rate  constant  ( k ),  they  simulated  cyclic  voltammetry  experi¬ 
ments  conducted  by  Uchida  et  al.  [20].  Four  distinct  peaks 
were  observed  with  two  each  on  the  anodic  and  cathodic 
scans.  Comparison  between  the  experimental  data  and  the 
model  predict  ions  yielded  values  of  2.2  x  1 0  9  cnr/s  and 
0.00019  cm5/2/(mol1/2  s)  for  and  k,  respectively.  Zhang 
et  al.  [19]  also  studied  the  effect  of  parameters  like  scan  rate 
on  the  predicted  cyclic  voltammogram.  They  observed  that 
when  the  scan  rate  is  increased,  the  height  of  all  the  peaks 
increased.  Further,  the  potentials  at  which  the  anodic  peaks 
occurred  shifted  to  more  positive  values  while  the  cathodic 
peaks  shifted  to  more  negative  values.  The  height  of  the  first 
peak  varied  linearly  with  the  square  root  of  scan  rate  and  the 
relationship  was  shown  to  agree  well  with  literature  values. 
The  other  peaks,  however,  could  not  be  related  linearly  to  the 
square  root  of  scan  rate  and  this  has  also  been  reported  by 
previous  researchers.  Zhang  et  al.  [19]  grouped  the  rate  of 
electrochemical  reaction  and  the  rate  of  lithium  diffusion 
into  a  single  dimensionless  parameter  b  given  by 


Keeping  the  value  of  D\  constant,  the  effect  of  varying  b  was 
studied.  Initially,  the  peak  current  increased  strongly  with  b 
but  after  a  point  the  peak  height  became  insensitive  to  b. 
Based  on  this  observation  the  authors  concluded  that  for  low 
values  of  b  both  kinetic  and  diffusion  limitations  were 
important  and,  therefore,  the  peak  current  increased  with 


b.  But  as  b  kept  increasing,  a  point  was  reached  when 
kinetics  became  too  fast  and  the  system  was  entirely  diffu¬ 
sion  limited.  As  a  result,  a  further  increase  in  b  did  not 
increase  the  peak  current.  The  authors  also  studied  the  effect 
of  the  anodic  transfer  coefficient  (aa)  on  the  simulated  cyclic 
voltammogram.  They  observed  that  for  high  values  of  b,  aa 
had  no  effect  on  the  cyclic  voltammogram  while  for  low 
values  of  b,  a  decrease  in  aa  favored  the  anodic  reaction 
more  and  made  the  values  of  peak  potentials  smaller.  This  is 
because  when  b  is  low,  kinetics  is  slow  and  consequently  a 
decrease  in  aa  has  an  effect  on  the  peak  potential.  On  the 
other  hand,  when  kinetics  is  fast,  i.e.  when  b  is  high,  change 
in  aa  does  not  produce  any  effect  on  the  voltammogram. 

The  models  developed  by  Doyle  et  al.  [1],  and  later 
extended  by  Pals  and  Newman  [4],  Arora  et  al.  [10,13] 
and  Botte  et  al.  [8],  were  mathematically  2D  with x  being  the 
dimension  across  the  cell  thickness  and  r  being  the  dimen¬ 
sion  along  the  radii  of  the  spherical  active  material  particles 
in  the  electrodes.  Since  r  does  not  represent  a  physically  new 
dimension,  it  is  often  referred  to  as  a  pseudo-dimension  and 
the  corresponding  models  as  pseudo-2D  models.  In  solving 
these  models  Doyle  et  al.  [1]  and  Arora  et  al.  [10,13]  used  a 
technique  called  Duhamel’s  super-position  to  collapse  the 
pseudo-dimension  r,  thereby  making  the  model  ID.  By  this 
method,  one  can  express  the  surface  concentration  on  the 
electrode  particle  directly  as  an  infinite  series  in  D\tlr 1 .  One 
disadvantage  of  this  method  is  that  in  order  to  obtain 
accurate  results  a  large  number  of  algebraic  calculations 
have  to  be  made  at  each  node  point  in  x  for  each  time 
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step.  This  calculation  can  be  reduced  significantly  by  using 
the  method  suggested  by  Subramanian  and  White  [21]  when 
solving  for  diffusion  in  the  solid  phase.  In  their  method,  a 
solution  for  local  concentration  (ci)  as  a  function  of  r  and  t  is 
assumed  to  be  of  the  form 

N 

ci  =  ffAi(t)r2‘  (41) 

;=o 

where  the  number  of  terms  (N)  depends  on  the  accuracy 
required.  The  coefficients  A,(f)  are  found  by  forcing  the  solu¬ 
tion  to  satisfy  the  boundary  conditions,  as  well  as  the  governing 
equation  at  conveniently  selected  points  in  r.  The  number 
of  parameters  A,(r)  to  be  determined  is  exactly  the  same  as 
N  and,  therefore,  the  corresponding  model  is  called  the 
^parameter  model.  Subramanian  and  White  [21]  conduct 
case  studies  for  three  different  values  of  N,  namely,  2,  3  and  4. 
As  a  result  of  comparison  the  exact  solution  with  the  two-, 
three-  and  four-parameter  models  for  constant  current  dis¬ 
charges,  the  authors  conclude  that  the  four-parameter  model 
predicts  solutions  with  good  accuracy  even  for  high  discharge 
rates  (say  10C).  The  two-  and  three -parameter  models,  how¬ 
ever,  make  significantly  erroneous  predictions  at  such  high 
rates  and  must  be  used  with  caution  even  at  moderate  rates. 
At  low  rates,  all  the  three  models  predict  results  that  match 
very  well  with  the  exact  solution.  This  method  was  used  in  a 
full  cell  model  by  Gomadam  et  al.  [7]  in  solving  for  lithium 
transport  in  the  active  material  solid  particles. 

1.7.  Conclusions 

A  review  of  the  modeling  efforts  of  researchers  in  the  field 
of  Li/Li-ion  batteries  at  the  University  of  South  Carolina 
is  presented.  The  review  included  models  that  describe 
the  behavior  of  single  electrode  particles,  single  electrodes 
and  full  cells  under  various  operating  conditions  like  cyclic 
voltammetry,  constant  current  discharge  and  impedance.  It 
was  concluded  that  the  models  presented  must  be  extended 
to  include  several  phenomena  like  concentration  dependen¬ 
cies  of  diffusion  coefficients,  state  of  charge  dependencies 
of  half  cell  entropies  and  salt  precipitation  at  high  discharge 
rates.  Efforts  are  currently  underway  at  the  University  of 
South  Carolina  to  understand  these  phenomena  so  as  to 
enable  our  models  to  predict  reality  closer. 

2.  Modeling  of  nickel  based  batteries 

2.1.  Background 

The  nickel  hydroxide  electrode  has  proved  to  be  one  of 
the  most  enduring  battery  chemistries,  with  widespread  use  in 
a  variety  of  applications  ranging  from  portable  electronics 
to  satellites.  Used  in  a  number  of  battery  systems,  including 
the  nickel-cadmium,  nickel-metal  hydride,  nickel-zinc 
and  nickel-hydrogen,  the  chemistry  in  the  nickel  hydroxide 


electrode  has  proved  sufficiently  pliable  to  suit  many  an 
application.  Present  interest  stems  from  use  in  space  applica¬ 
tions,  including  in  the  Hubble  space  telescope  (Ni-H2)  [22] 
and  for  use  in  electric  cars  (Ni-MH)  [23].  The  original  patent 
for  the  nickel-cadmium  cell  can  be  traced  back  to  the  last 
decade  of  the  19th  century  and  is  attributed  to  a  Swedish 
scientist,  lunger  [24]  (cited  by  Salkind  [25]).  Its  first  commer¬ 
cial  use  in  a  large  scale  was  reported  in  1917,  when  it  was  used 
to  light  subway  trains  in  Paris  [26].  Since  then,  the  advantages 
of  the  battery  chemistry  (e.g.  cycle  life,  shelf  life  etc.)  has 
ensured  its  dominance  of  the  secondary  battery  market  [27]. 

However,  a  century  of  research  and  development  into 
this  chemistry  has  not  resulted  in  a  fail-safe  method  of 
fabricating  cells.  Part  of  the  problem  stems  from  the  com¬ 
plex  nature  of  the  nickel  hydroxide  redox  behavior  and  the 
close  relationship  between  the  structure  of  the  material 
and  its  electrochemical  characteristics  [28].  In  addition,  the 
miniaturization  of  most  modern  devices  has  thrown  up  new 
challenges  in  the  field  of  design  of  batteries.  Therefore,  one  is 
confronted  with  not  only  understanding  the  various  physio- 
chemical  processes  that  occur  in  battery  systems,  but  also  in 
designing  them  to  suit  our  ever  increasing  demands.  Math¬ 
ematical  models  derived  from  continuum  mechanics  provide 
an  effective  tool  for  achieving  both  these  ends.  With  the 
development  of  modern  computers  with  enhanced  speed, 
enormous  opportunities  are  available  for  the  development 
of  robust  theoretical  descriptions  of  complex  systems.  In  this 
paper,  we  detail  the  efforts  that  have  been  undertaken  to 
model  the  nickel  hydroxide  active  material  over  the  last  two 
decades.  The  purpose  of  the  review  is  two-fold  namely:  (i)  to 
give  the  reader  an  understanding  of  the  theoretical  basis  for 
the  various  electrochemical  signatures  observed  in  the  nickel 
electrode  and  (ii)  to  detail  the  efforts  in  developing  design 
tools  that  relate  various  electrode  properties  (e.g.  thickness, 
porosity)  to  the  performance. 

2.2.  System  description 

Fig.  1  schematizes  a  cross-section  of  a  cell  constructed 
with  the  nickel  positive  electrode.  The  negative  electrode 
used  in  the  cell  could  be  a  porous  electrode  (Cd,  MH)  or  a 
planar  electrode  (H2).  The  nickel  hydroxide  active  material 
is  supported  on  a  conductive  matrix  with  the  active  material 
either  pasted  or  impregnated  into  the  pores  of  the  sinter. 
There  is  considerable  void  volume  in  the  electrode,  which  is 
flooded  with  electrolyte,  typically  concentrated  KOH  in  the 
concentration  range  of  26-33%  [29].  The  electrodes  repre¬ 
sented  in  the  figure  can  be  thought  of  as  half  the  actual 
electrode  thickness,  with  equivalent  processes  occurring  in 
the  other  half.  Symmetry  enables  us  to  neglect  the  processes 
in  the  other  half.  The  main  reaction  that  occurs  at  the  nickel 
electrode  can  be  ideally  represented  as  [30] 

charge 

Ni(OH)2  +  OH  NiOOH  +  H20  +  e", 

discharge 

Eq  =  0.311  V  versus  Ag/AgCl  (42) 
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A  proton  and  electron  are  released  at  the  nickel  site  within 
the  solid  active  material.  The  electron  travels  to  the  con¬ 
ducting  substrate  and  eventually  out  to  the  external  circuit. 
The  proton  diffuses  through  the  active  material  to  the  solid/ 
electrolyte  interface  where  it  combines  with  the  hydroxyl 
ion  to  form  water.  On  discharge,  these  processes  are 
reversed. 

On  both  charge  and  discharge,  oxygen  evolution  occurs 
simultaneously  with  reaction  1,  as  the  equilibrium  potential 
of  this  reaction  is  more  negative. 

40H~  ->  02  +  2H20  +  4e  , 

E0  =  0.22  V  versus  Ag/AgCl  (43) 

In  a  sealed  cell,  the  oxygen  evolved  in  the  positive  electrode 
migrates  to  the  negative  electrode  and  is  reduced  (reverse  of 
reaction  (43)),  thus  ensuring  that  no  excessive  pressure 
builds  up  in  the  cell.  The  oxygen  transfer  in  the  cell  is 
facilitated  by  using  an  oxygen  permeable  membrane  and 
operating  the  cell  under  starved  conditions,  i.e.  with  limited 
amount  of  electrolyte. 

Again,  reaction  (42)  is  a  simplistic  representation  of  the 
redox  processes  that  occur  in  the  electrode,  as  the  active 
material  is  known  to  consist  of  different  phases.  Bode  et  al. 
[30]  summarized  the  complex  chemical  and  electrochemical 
processes  as  follows. 


a  -  Ni(OH)2  y  -  NiOOH 

*  discharge 


Dehydration 


Overcharge  (44) 


P  -  Ni(OH)2  .  p- NiOOH 

L  discharge 


Note  that  these  reactions  are  not  written  as  balanced  reac¬ 
tions,  but  rather  to  illustrate  the  relationship  between  the 
various  phases.  The  hydrated,  non-close-packed  phase, 
termed  a-Ni(OH)2,  is  unstable  and  dehydrates  to  the  anhy¬ 
drous  close-packed  phase,  termed  (3-Ni(OH)2,  in  concen¬ 
trated  alkali  solution.  Oxidation  of  |3-Ni(OH)2  results  in  the 
formation  of  P-NiOOH  while  oxidation  of  a-Ni(OH)2  results 
in  v-NiOOH.  The  P-phase  converts  to  the  y-phase  on  over¬ 
charge,  as  Ni(III)  is  converted  to  Ni(IV).  Loyselle  et  al.  [62] 
used  a  point  defect  approach  to  modify  this  reaction  mecha- 
mism. 


2.3.  Cell  models 

In  order  to  predict  the  performance  of  a  nickel  battery  for 
a  wide  range  of  designs  and  operating  conditions,  a  number 
of  processes  must  be  accounted  for  in  terms  of  measurable 
parameters:  (i)  thermodynamics  and  kinetics  of  all  redox 
reactions;  (ii)  transport  of  ions  through  all  ionically 
conducting  phases;  (iii)  the  potential  variations  in  both 
ionically  and  electronically  conducting  phases.  The  high 
concentration  of  OH-  in  the  cell  (31%)  necessitates  the 


use  of  concentrated  solution  theory  [31-35],  where  the 
interactions  of  the  ions  with  the  solvent  also  need  to  be 
considered.  The  complications  introduced  by  the  presence 
of  the  porous  matrix  are  handled  using  porous  electrode 
theory  [34,36-38], 

Macroscopic,  porous  electrode  models  have  been  devel¬ 
oped  for  the  nickel  batteries  [39-44]  in  order  to  understand 
how  the  complex  interactions  that  occur  among  the  anode, 
cathode  and  separator  affect  battery  performance.  For  exam¬ 
ple,  Micka  and  Rousar  [39,40]  developed  a  ID  model,  but 
they  noticed  that  on  discharge  the  voltage  curves  could  not 
be  fitted  adequately  using  normal  Nernstian  thermody¬ 
namics,  especially  at  high  states  of  charge.  Fan  and  White 
[43]  improved  the  predictions  of  the  Ni-Cd  discharge  curves 
by  using  activity  coefficient  corrections  to  the  Nernst  equa¬ 
tion,  as  suggested  by  Barnard  et  al.  [45].  However,  their 
model  was  unable  to  adequately  predict  the  percent  utiliza¬ 
tion  of  the  active  material  on  discharge  as  a  function  of 
current.  This  limitation  resulted  from  assuming  that  the 
nickel  hydroxide  active  material  charged  and  discharged 
uniformly. 

De  Vidts  and  White  [46]  modeled  a  sealed  Ni-Cd  cell  that 
includes  proton  diffusion  and  ohmic  drop  through  the  active 
material  in  the  nickel  electrode.  The  model  was  used  to 
calculate  the  sensitivity  coefficients  for  various  parameters 
in  the  model.  These  calculations  showed  that  the  discharge 
voltage  of  the  cell  is  affected  mostly  by  the  kinetics  of  the 
nickel  reaction.  Towards  the  end  of  discharge,  proton  diffu¬ 
sion  was  also  found  to  become  important,  because  the  proton 
diffusion  process  affects  the  active  material  utilization  sig¬ 
nificantly.  During  charge,  the  cell  voltage  is  mainly  affected 
by  the  kinetics  of  the  nickel  reaction  until  the  oxygen 
evolution  reaction  begins,  after  which  time  the  kinetics  of 
the  oxygen  evolution  has  the  largest  effect.  The  oxygen 
evolution  reaction  is  also  the  most  influencing  factor  on  the 
actual  charge  uptake  of  the  cell  by  the  end  of  a  charge 
operation  (charge  efficiency).  Compared  with  the  rates  of 
reaction  and  proton  diffusion,  the  ohmic  drop  in  the  active 
material  of  the  nickel  electrode  and  the  mass  transport  and 
ohmic  drop  in  the  electrolyte  have  negligible  effect  on  the 
behavior  of  the  cell  studied  here. 

De  Vidts  et  al.  [47]  also  developed  a  model  for  predicting 
the  discharge  behavior  of  metal  hydride  electrode.  The 
model  was  used  to  study  the  effect  of  various  parameters 
on  predicted  discharge  curves.  The  simulations  obtained 
using  the  model  show  the  expected  decrease  of  charge 
utilization  as  the  rate  of  discharge  is  increased.  Increasing 
the  particle  size  of  the  alloy  and  decreasing  the  diffusion 
coefficient  of  the  hydrogen  atoms  in  the  hydride  showed  a 
similar  effect  on  the  discharge  curves.  The  model  simula¬ 
tions  also  show  the  critical  role  that  the  kinetic  and  transport 
parameters  play  in  determining  the  overall  shape  of  the 
predicted  discharge  curves  for  a  metal-hydride  electrode. 

Mao  and  White  [48]  used  a  mathematical  model  to 
characterize  the  self-discharge  of  a  NiOOH  electrode  in  a 
hydrogen  environment.  The  model  included  diffusion  of 
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dissolved  hydrogen  in  an  electrolyte  film,  which  covers  a 
flooded  electrode,  electrochemical  oxidation  of  hydrogen, 
reduction  of  NiOOH,  and  changes  of  surface  area  and 
porosity  of  the  electrode  during  the  self-discharge  process. 
Although  the  self-discharge  process  is  complicated,  the 
predictions  of  the  model  are  consistent  with  experimental 
results  reported  in  the  literature,  which  include  linear  rela¬ 
tionships  between  the  logarithm  of  hydrogen  pressure  and 
time  and  between  the  logarithm  of  the  capacity  remaining 
and  time.  The  model  predictions  indicate  that  hydrogen 
oxidation  takes  place  predominantly  near  the  front  side  of 
the  electrode,  but  the  reduction  of  NiOOH  takes  place 
uniformly  throughout  the  electrode. 

Mao  et  al.  [49]  improved  the  percent  utilization  predic¬ 
tions  of  Barnard  et  al.  [45]  by  adding  the  diffusional 
resistance  of  protons  to  their  pseudo  2D  model.  However, 
this  model  is  unable  to  predict  the  change  in  performance  of 
the  battery  on  cycling  because  the  redox  reaction  in  the 
active  material  is  assumed  to  be  the  overly  simplistic  reac¬ 
tion  (1).  Adding  a  more  sophisticated  reaction  mechanism  to 
the  model  is  required  in  order  to  predict  performance  during 
cycling. 

During  the  discharge  of  nickel  hydroxide  electrodes,  two 
voltage  plateaus  have  been  frequently  observed  [29,50]. 
While  the  first  plateau  occurs  due  to  the  nickel  oxidation/ 
reduction  reaction  via  reaction  (42),  the  second  plateau 
appears  approximately  400  mV  cathodic  of  it  and  has  been 
reported  to  have  as  much  as  50%  of  the  capacity  of  the  first 
plateau  [29].  Mao  et  al.  [49]  modeled  the  second  discharge 
plateau  in  a  nickel-hydrogen  cell  using  a  comprehensive  cell 
model,  which  incorporates  proton  diffusion  and  porous 
electrode  theory  using  a  pseudo  2D  approach.  They  con¬ 
sidered  the  conductivity  of  the  active  material  to  be  a 
continuous  function  until  a  state  of  change  of  0.1.  At  this 
point,  the  authors  assumed  a  discontinuity  in  the  conduc¬ 
tivity  at  a  constant  low  value.  This  can  be  interpreted  as  a 
model  that  incorporates  the  first  explanation  of  the  second 
plateau  given  where  the  conductivity  of  the  active  material 
decreases  steadily  during  discharge  between  the  state  of 
charge  of  1  and  0.1.  At  this  low  state  of  charge  the  conduc¬ 
tion  in  the  film  shifts  from  electron  conducting  to  hole 
conducting  and  the  hole  conductivity  is  a  constant  at  around 
1.2  x  10  8/(£2  cm).  This  low  conductivity  results  in  a 
greater  ohmic  drop  and  hence  the  potential  where  the 
reaction  occurs  is  lower. 

In  contrast  to  the  models  that  assume  the  conductivity  to 
be  the  cause  for  the  second  discharge  plateau,  the  model  by 
De  Vidts  et  al.  [51]  was  successful  in  predicting  the  second 
discharge  plateau  for  a  Ni-H2  cell  by  considering  the  oxygen 
reduction  reaction  (reverse  of  reaction  (43)).  Although,  the 
oxygen  reduction  reaction  has  been  included  in  previous 
models  for  the  Ni-Cd  battery  developed  by  the  authors 
[42-44,52]  they  do  not  predict  the  second  plateau.  The 
authors  explain  this  by  suggesting  that  the  oxygen  generated 
in  the  positive  electrode  in  a  Ni-Cd  cell  is  transported  to  the 
negative  electrode  easily  as  the  mass  transfer  resistance  is 


low.  At  the  negative  electrode,  the  oxygen  is  reduced.  Hence, 
not  much  dissolved  oxygen  is  present  in  these  cells.  How¬ 
ever,  in  a  Ni-H2  cell,  the  model  predicts  that  there  is  ample 
oxygen  present  to  sustain  the  reaction,  which  appears  to  be 
due  to  the  gas  phase  present.  De  Vidts  et  al.  [51]  simulated 
the  effect  of  oxygen  pressure  on  the  second  discharge 
plateau  and  predict  a  change  in  potential  to  less  negative 
values  on  increasing  the  pressure.  In  addition,  the  models 
predict  that  an  increase  the  pressure  would  result  in  a  bigger 
second  plateau,  as  there  is  more  oxygen  in  the  cell. 

The  model  developed  by  De  Vidts  et  al.  [51]  was  extended 
by  Wu  et  al.  [53]  to  include  phase  reactions  in  the  nickel 
active  material  of  a  Ni-H2  cell.  The  phase  reactions  occurred 
in  two  parallel  redox  paths  with  different  electron  transfer 
numbers:  (i  -  N  i  ( O  H)  2/  [1  -  N  i O O H  with  1.0  electron  transfer 
and  a-Ni(OH)2/y-NiOOH  with  1.67  electron  transfer  and 
two  connection  reactions  between  them:  a-Ni(OH)2/p- 
Ni(OH)2  and  p - N i O O H /y - N i O OH.  Important  mechanisms 
inside  a  Ni-H2  cell,  such  as  mass  balances  of  active  species, 
kinetics  of  electrochemical  reactions,  and  the  energy  bal¬ 
ance  of  the  whole  cell,  etc.  have  been  considered  in  the 
model.  The  model  predictions  under  different  conditions 
were  presented  and  compared  to  experimental  data.  The 
predicted  cell  voltage  and  cell  pressure  agreed  well  with  the 
experimental  data  as  shown  in  Fig.  8.  Good  agreement 
between  experiment  and  model  was  also  obtained  with 
the  temperature-time  curves  (not  shown  here).  Based  on 
the  simulation  results,  it  was  concluded  that  nickel  phase 
reactions  have  significant  influences  on  the  behavior  of  a  Ni- 
H2  cell.  Further,  some  observed  phenomena  like  the  capacity 
variation  at  different  temperatures,  KOH  concentration 
change  between  charge  and  discharge  and  cell  potential 
rollover  during  overcharge  were  explained  with  the  model. 

Wang  et  al.  [54]  developed  a  micro-macroscopic  coupled 
model  of  electrochemical  kinetics  and  transport  phenomena 
occurring  in  batteries  and  fuel  cells  using  the  technique  of 
volume  averaging.  The  model  accounted  for  the  effects  of 
the  micro-scale  and  interfacial  non-equilibrium  processes  on 
the  macroscopic  species  and  charge  transfer.  It  was  shown 
that  solid-state  diffusion,  ohmic  resistance  in  resistive  active 
materials  and  electrolytes,  and  interface  passivation  can  all 
be  incorporated  in  the  present  model  in  a  consistent  manner. 
This  ultimately  allowed  for  the  integration  of  detailed 
chemistry  and  morphology  of  the  electrode/electrolyte  inter¬ 
face  into  a  battery  model.  Gu  et  al.  [55]  applied  this  model 
to  predict  discharge  and  charge  behaviors  of  Ni-Cd  and 
Ni-MH  cells.  The  model  predictions  were  validated  against 
the  previous  results  available  in  the  literature  for  a  Ni-Cd 
cell  and  a  single  MH  electrode  with  good  agreement.  As 
compared  with  the  pseudo  2D  numerical  model,  the  present 
model  offers  a  more  efficient  approach  to  the  modeling 
of  intercalative  batteries  and  can  be  more  easily  applied  to 
simulate  complex  cycles  of  discharge,  rest  and  recharge  as 
involved  in  the  electric  vehicle  application.  The  results  of  the 
simulation  showed  that  the  Ni-MH  cells  behave  somewhat 
differently  from  Ni-Cd  cells.  In  particular,  it  was  observed 


PM.  Gomadam  et  al.  /  Journal  of  Power  Sources  110  (2002)  267-284 


281 


that  both  hydrogen  diffusion  in  MH  particles  and  proton 
diffusion  across  the  nickel  active  material  can  be  dominant 
mechanisms  limiting  cell  performance,  depending  on  the 
ratio  of  charge  capacities  of  the  two  electrodes. 

2.4.  Modeling  of  the  Ni(OH)2  active  material 

As  the  aforementioned  examples  illustrate,  the  limitation 
of  previous  nickel  battery  models  has  been  due  to  the  limited 
understanding  of  how  the  nickel  hydroxide  active  material 
behaves.  As  improvements  have  been  made  in  predicting  the 
behavior  of  the  active  material,  the  predictions  of  the  battery 
models  have  improved.  What  follows  is  a  review  of  models 
that  have  been  developed  to  predict  various  phenomena  that 
occur  in  the  nickel  hydroxide  active  material.  This  review 
will  include  models  that  describe  corrections  to  the  Nernst 
equation,  proton  diffusion,  and  mechanisms  of  the  redox 
process,  including  hysteresis.  Ultimately,  the  information 
that  is  gained  from  models  of  the  active  material  needs  to  be 
incorporated  into  models  that  couple  variations  in  the  micro¬ 
scopic  active  material  with  variations  in  the  macroscopic 
porous  electrodes. 

2.4.1.  Corrections  to  the  Nernst  equation 

The  equilibrium  potential  of  the  nickel  electrode  depends 
on  the  state  of  charge  of  the  electrode  as  governed  by  the 
Nernst  equation.  However,  investigations  have  shown  that 
for  a  greater  part  of  the  oxidation  process  the  potential  is 
independent  of  the  degree  of  oxidation  [29],  a  feature  that 
cannot  be  described  by  the  Nernst  equation.  In  particular,  as 
pointed  out  by  Paxton  and  Newman  [56],  two  main  features 
define  the  discharge  curves,  namely:  (i)  during  most  of  the 
discharge,  the  curve  appears  Nernstian  and  (ii)  at  low  states 
of  charge  the  potential  change  is  steeper  than  that  predicted 
by  Nernstian  thermodynamics.  This  deviation  from  Nerns¬ 
tian  thermodynamics  also  occurs  at  high  states  of  charge 
[57].  Barnard  et  al.  [45]  studied  the  reversible  potential  of 
the  reaction  by  re-deriving  the  Nernst  equation  using  activ¬ 
ities  in  place  of  concentrations.  They  studied  the  effect  of  the 
one-parameter  Margules  and  two-parameter  Van  Laar  activ¬ 
ity  coefficient  models  on  the  reversible  potential.  With  these 
corrections  the  equilibrium  potential  takes  the  form 

E  =  E°  +  —  Inf—* 9— )  +  —  [— 2A  +  4A(1  -  9)  -  2B 
nF  \l-0j  2Fl  1  ; 

+  6£(1  -B)  -3£(1  -6)2}  (46) 

with  the  expression  simplifying  to  the  Nernst  equation  for 
A  =  B  =  0. 

The  authors  demonstrated  that  changes  in  the  parameters 
(A  and  B )  can  be  used  to  predict  the  discharge  profile  of  the 
nickel  hydroxide  electrode  and  suggested  that  the  potential 
of  the  electrode  was  held  constant  by  “pairs  of  co-existing 
solid  solutions”.  Jain  et  al.  [57]  recently  extended  Barnard’s 
model  by  considering  more  models  for  the  activity  coeffi¬ 
cient.  They  compared  their  model  to  experimental  data  on 


thin  films  of  nickel  hydroxide  with  12%  cobalt  additive  and 
extracted  the  activity  coefficient  parameters  and  the  equili¬ 
brium  potential  as  a  function  of  temperature.  The  parameters 
extracted  were  for  a  cathodically  deposited  film,  which  has 
been  cycled  a  few  times  and  would  represent  the  a-y-phase. 
Films  that  have  a  different  structure  would  require  a  different 
set  of  parameters  to  represent  them.  The  authors  also 
suggested  that  the  insertion  of  protons  into  the  lattice 
occurred  in  an  ordered,  non-random  fashion. 

Fig.  7  shows  the  equilibrium  potential  curves  obtained  by 
plotting  the  functionality  in  Eq.  (46).  The  voltage  has  been 
referenced  to  the  voltage  of  the  nickel  hydroxide  electrode  at 
50%  state  of  charge.  Fig.  7  shows  the  voltages  as  described 
by:  (a)  the  Nernst  equation  (A  =  B  =  0  in  Eq.  (46));  (b)  the 
one-parameter  model  (A  =  1.198,  B  =  0);  (c)  the  two-para- 
meter  model  (A  =  4.44,  B  =  —4.576),  the  value  backed  out 
by  Jain  et  al.  [57]  for  nickel  hydroxide  films  with  12%  cobalt 
additive,  at  27  °C.  Increasing  the  intercalation  constant 
flattens  the  voltage  profile  through  most  of  the  discharge. 
In  addition,  the  sharp  drop  in  potential  at  low  and  high  states 
of  charge  is  attenuated.  Both  these  features  have  been  seen 
experimentally  in  the  nickel  hydroxide  electrode  [57].  As 
expected,  the  curve  follows  the  Nernst  equation  through 
most  states  of  charge. 


2.4.2.  Proton  diffusion 

The  models  that  predict  the  constant  current  discharge  of 
the  cell  show  that  the  value  used  for  the  diffusion  coefficient 
of  protons  is  critical  in  predicting  the  utilization  at  high 
discharge  rates.  Weidner  and  Timmerman  [58]  were  the  first 
to  show  the  importance  of  the  parameter  in  terms  of  its  effect 
on  the  utilization  of  the  active  material.  Later,  Mao  et  al.  [49] 
showed  that  while  using  a  value  of  1  x  1CV10  cm2/s  ensured 
complete  utilization  for  a  10C  discharge,  using  1  x  10  12 
cnf/s  results  in  an  80%  loss  in  utilization.  Unfortunately, 
early  attempts  [44]  to  evaluate  the  proton  diffusion  coefficient 
resulted  in  values  that  varied  by  three  orders  of  magnitude 
from  1  x  10  9  to  1x10  12  cm2/s.  Therefore,  Motupally 
et  al.  [59,60]  used  electrochemical  impedance  spectroscopy, 
with  the  aid  of  a  first-principles  model,  to  measure  the 
diffusion  coefficient  of  protons  in  a  thin  film  of  nickel 
hydroxide  with  12%  cobalt  additive  deposited  on  a  gold 
substrate.  They  evaluate  the  slope  of  the  impedance  from 
their  equations  and  equate  it  to  the  experimental  data  in  the 
transition  region  (transition  from  45°  to  slope  of  oo)  to 
estimate  the  diffusion  coefficient  as  a  function  of  state  of 
charge.  The  authors  found  that  the  diffusion  coefficient  data 
changed  as  a  function  of  state  of  charge  and  followed  the  trend 
given  by 


Dh+  —  Z)n;ooh 


g+(l_g)fDNi(OH)A 

V  ^NiOOH  J 


1/2' 


(47) 


with  £>Ni(0H)2  =  6.4  x  HU11  cm2/s  and  ONi00H  =  3.4x 
1 0  K  cm2/s.  Eq.  (47)  is  a  mixing  rule  in  terms  of  the  root 
mean  square  displacement  of  the  diffusing  species  in  a  solid 
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solution  comprising  of  a  homogeneous  mixture  of  two 
phases  [60]. 

Motupally  et  al.  [61]  also  extended  the  constant-current 
discharge  model  by  Weidner  and  Timmerman  [58]  to 
include  the  state  of  charge  dependence  of  the  diffusion 
coefficient.  They  compared  the  two  models  and  found  that 
while  on  charge,  the  constant  diffusion  model  compares 
adequately  with  the  exact  solution,  there  is  a  significant 
deviation  on  discharge.  Using  the  model,  the  authors  predict 
the  utilization  of  the  active  material  as  a  function  of  the 
charge  and  discharge  currents.  They  then  verify  the  validity 
of  the  proton  diffusion  functionality  with  state  of  charge 
(Eq.  (47))  by  conducting  experiments  on  thin  films  of  nickel 
hydroxide.  The  utilization  of  the  active  material  was  found 
to  be  close  to  100%  irrespective  of  the  C-rate  on  charge, 
while  there  is  a  significant  decrease  in  utilization  on  dis¬ 
charge.  The  excellent  agreement  between  the  model  and  the 
experiments  gives  credence  to  the  diffusion  coefficients 
estimated  by  Motupally  et  al.  [60].  It  also  validates  the 
assumption  that  the  rate  of  the  nickel  hydroxide  electrode  is 
governed  by  transport  limitations,  as  the  model  by  Motu¬ 
pally  et  al.  [61]  does  not  include  the  kinetics  of  the  nickel 
reaction.  The  authors  also  attempted  to  estimate  a  constant 
diffusion  coefficient  that  would  adequately  fit  the  data  at  all 
C-rates.  However,  the  constant  diffusion  coefficients  used 
are  highly  inaccurate  in  predicting  the  utilization  on  dis¬ 
charge  and  under-predict  the  utilization  on  charge. 

2.4.3.  Reaction  mechanism 

Loyselle  et  al.  [62]  used  a  point  defect  approach  to 
describe  the  active  material  and  proposed  that  the  difference 


between  the  two  phases  was  the  level  of  defects.  They 
described  the  various  phases  in  the  active  material  using 
the  concentration  of  protons,  alkali  cations,  Ni  vacancies, 
and  oxygen  present  in  the  lattice,  and  showed  the  applic¬ 
ability  of  the  defect  model  in  describing  some  of  the 
observed  phenomenon  (e.g.  existence  of  a  maximum  oxida¬ 
tion  state  of  3.67).  With  this  information  the  authors  mod¬ 
ified  the  Bode  diagram  to  incorporate  these  features  (for 
more  details  on  the  non-stoichiometric  structural  model  and 
the  modified  Bode  diagram  refer  to  (44)).  Using  this  mod¬ 
ified  Bode  diagram,  Srinivasan  et  al.  [63]  developed  a 
procedure  that  relates  the  defect  parameters  to  the  mass 
change  and  capacity.  The  procedure  was  used  to  understand 
four  frequently  observed  phenomena  in  the  nickel  electrode, 
namely  (i)  1.67e_  transfer  in  the  first  charge  while  subse¬ 
quent  charges  lead  to  le_  transfer;  (ii)  le_  transfers  for  the 
first  discharge;  (iii)  a  steady  decrease  in  capacity  on  both 
charge  and  discharge  on  cycling;  (iv)  an  increase  in  the  total 
mass  of  the  film  on  cycling. 

Hysteresis  is  another  phenomenon  that  greatly  affects  the 
performance  of  nickel  batteries.  The  nickel  hydroxide  elec¬ 
trode  is  known  to  exhibit  a  stable  hysteresis  loop,  with  the 
potential  on  charge  being  higher  than  that  on  discharge  at 
every  state  of  charge.  However,  very  little  beyond  this 
observation  has  been  quantified.  Recently,  Srinivasan, 
et  al.  [64]  showed  that  this  hysteresis  loop  created  during 
a  complete  charge  and  discharge  (i.e.  boundary  curves)  is 
not  sufficient  to  define  the  state  of  the  system.  Rather, 
internal  loops  within  the  boundary  curves  (i.e.  scanning 
curves)  can  be  generated  that  access  potentials  between 
the  boundary  curves.  The  potential  obtained  at  any  state 
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of  charge,  as  well  as  how  the  material  charges  and  discharges 
from  that  point,  depends  on  the  cycling  history  of  the 
material.  The  implication  of  this  phenomenon  is  that  the 
potential  of  nickel-based  batteries  cannot  be  used  as  an 
indication  of  the  state  of  charge  of  the  cell.  Analysis  of 
the  boundary  and  scanning  curves  suggest  that  the  electrode 
consists  of  a  number  of  individual  units  or  domains,  each  of 
which  exhibits  two  or  more  metastable  states.  They  dis¬ 
cussed  the  cycling  behavior  of  the  nickel  hydroxide  elec¬ 
trode  within  the  context  of  previously  developed  theoretical 
arguments  regarding  domain  theory.  Although  the  specific 
cause  for  the  metastability  in  each  domain  is  not  understood, 
considerable  insights  are  provided  into  the  history-depen- 
dent  behavior  of  the  nickel  hydroxide  electrode.  Finally, 
they  developed  an  empirical  procedure  to  predict  the  scan¬ 
ning  curves  based  on  the  boundary  curves. 


2.5.  Conclusions 

A  review  of  the  continuum  models  developed  at  the 
University  of  South  Carolina  to  describe  the  nickel  battery 
systems  was  presented.  The  review  also  includes  models  that 
describe  detailed  behavior  of  the  nickel  hydroxide  active 
material.  The  methodology  used  in  developing  the  models, 
and  their  relative  strengths/weaknesses  were  discussed. 
Future  improvements  to  these  models  will  come  through 
a  better  understanding  of  how  the  active  material  cycles. 
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